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Matrix tablets are drug delivery devices designed to release a drug in a controlled manner over 
an extended period of time. We develop a cellular automaton (CA) model for the dissolution and 
release of a water-soluble drug and excipient from a matrix tablet of water-insoluble polymer. Cells 
of the CA are occupied by drug, excipient, water or polymer and the CA updating rules simulate 
the dissolution of drug and excipient and the subsequent diffusion of the dissolved substances. In 
addition we simulate the possible fracture of brittle drug and excipient powders during the tablet 
compression and the melting of the polymer during a possible thermal curing process. Different 
stirring mechanisms that facilitate the transport of dissolved drug in the fluid in which the tablet is 
immersed are modeled in the water cells adjacent to the boundary of the tablet. We find that our 
simulations can reproduce experimental drug release profiles. Our simulation tool can be used to 
streamline the formulation and production of sustained release tablets. 
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The design of controlled release systems has been an active area of research in pharmaceuti- 
cal science and industry for decades. The most common systems include coated systems, matrix 
tablets, eroding tablets and oral osmotic therapeutic systems [TH3]- The present work focuses on 
the development of a mathematical model and simulation tool that describes the sustained drug 
release observed in matrix tablets. Matrix tablets are devices that deliver a drug in a controlled 
manner over an extended period of time. A preferred manufacturing method is to mix the drug with 
pharmaceutically inactive excipient and polymer powders. This powder mixture is then compressed 
in a die at high pressures (e.g. 70 — 200 MPa) and may be cured for 8-24 hours at 40 — 70 °C. Until 
now, pharmaceutical scientists were restricted to fabrication of experimental tablets to understand 
the influence of parameters such as the powder composition, the compaction pressure and the curing 
temperature and duration. As this is expensive and time consuming, the need for mathematical 
modeling and simulation tools becomes evident, see [1HZ] for earlier works. The aim of this work 
is to provide a tool that allows pharmaceutical scientists to mimic the release processes in matrix 
tablets. 

In our earlier work [7] we proposed two mathematical models for the release of a water-soluble 
drug from a polymer/excipient matrix tablet. The first model used a biased random walk on the 
contact graph of a random sphere packing. The second model was based on a system of reaction- 
diffusion partial differential equations for the concentrations of dissolved and undissolved drug and 
excipient, respectively. While the first discrete model predicted partial release of drug from the 
tablet in agreement with experimental observations [7J, the second continuous model proved better 
at capturing a change from convex to concave in several experimental release profiles. With the 
cellular automaton (CA) model presented here, we explicitly model the initial wetting process of the 
tablet when water passes through pores into the interior of the tablet, followed by the dissolution 
of the drug and the excipient. We also incorporate the fracture of the brittle drug and excipient 
powders during the compaction process and the partial melting of the polymer during an optional 
thermal treatment of the compressed tablet. This results in the formation of larger connected 
regions of polymer that may entrap drug and prevent it from leaving the tablet. Finally, we also 
allow to simulate different stirring protocols that aide the transport of dissolved and released drug 
away from the tablet [8]. We find that our model has the capacity to reproduce quantitatively 
experimental release curves of tablets formulated from different powder mixtures and subjected to 
different heating protocols. 

Cellular automata have shown their usefulness in many sciences since their introduction by Ulam 
and von Neumann about 70 years ago. Let us mention here only some applications to traffic 
networks [H] , neural networks (TO] , tumor growth [TT] and statistical mechanics fT2] . They allow for 
greater detailed modeling of internal processes in "cells" that is generally not possible in coarse- 
grained partial differential equation models. The recent advances in computational technology have 
made earlier restrictions to small numbers of cells less stringent, if not obsolete. 

The remainder of this paper is organized as follows. We describe in Section [iTj the cellular 
automaton model in detail, paying special attention to the initialization process. Section |TO| contains 
several simulation results, where we explore the influences of the respective parameters. In Section 
we discuss our results and indicate future research. Appendix [A] lists parameter values used in the 
simulations and instructions how to download and use our simulation tool. 
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II. THE CELLULAR AUTOMATON MODEL 

A. General description and program flow 

The cellular automaton recreates the three-dimensional cylindrical geometry of the matrix tablet. 
Each particle in the tablet is represented by a vertex in a cubic lattice that fills out a cylindrical 
domain. The diameter of the cylinder and height as a fraction of the diameter are specified by the 
user. A vertex interacts with its 6 closest von Neumann neighbors by a set of rules designed to 
imitate the actual physics of the release process, see Figure [TJ 




FIG. 1: (Left) The tablet is represented by a cubic lattice in a cylindrical domain. Each grid 
element, indicated here by a colored cube, can be in one of five states: drug (purple), polymer 
(red), excipient (green), empty (black) and water (blue). This tablet is surrounded by a layer 
of water cells (not shown here). (Right) Each element interacts with its 6 closest von Neumann 

neighbors. 

Every cell can be in one of five states, namely drug (D), polymer (P), excipient (X), empty (E) 
and water (W). In the initialization phase the cells of the lattice are randomly assigned a state so 
that the distribution of cells of type D, P and X are those prescribed by the composition of the 
powder mixture, with the possible formation of a polymer shell under thermal treatment. Empty 
cells are created during the heating and compression phase. Cells of type W surround the tablet at 
all times. The wet cells are further characterized by two quantities, the concentration of drug and of 
excipient, respectively. These quantities are real numbers between and 1, where a concentration 
of 1 corresponds to a saturated solution. Naturally, the concentrations change as water travels 
into the tablet and solid particles dissolve. To avoid clumsy triple index notation, we simply write 
CD(x,t) and cx(x,t) for the concentration of drug and excipient, respectively, in the cell located 
at vertex x of the lattice at time t. Throughout the paper, let N(x) denote the von Neumann 
neighborhood of vertex x. 

After initialization of the tablet according to the user's specification, the tablet's state is repeat- 
edly updated until either a maximum number of iterations have been executed or the amount of 
released drug ceases to change significantly. Each iteration of the update rules advances the sim- 
ulation clock by a certain fixed amount of time which depends on the model parameters. Finally, 
we also model the transport of dissolved drug in the boundary layer surrounding the tablet. The 
details of the simulation are described in the following sections. 
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FIG. 2: Schematic representation of the flow. During each update of the tablet, a new tablet state 
is created from the existing one which is afterwards cleared to free memory. 



B. Initialization 



The composition of the tablet is determined by the fractions of drug, polymer and excipient 
in the powder, which we denote here by P]j,Pp and Px, respectively. Thermal treatments and 
compaction studies on Eudragit-RLPO matrix tablets revealed that the particles of the powder 
mixture move and diffuse during the heating process |13l 114) , coalescing the polymer particles and 
forming matrices with decreased porosity and increased tortuosity. Porosity and tortuosity are 
important factors that affect the dissolution and release rate, as well as the possible formation of 
a thin polymer film on the surface of the tablet [13] . The release kinetics from laboratory studies 
show that, as the amount of polymer is increased, drug release from heated tablets approaches 
"first-order kinetics", i.e. nearly linear release with respect to time, see [7J Figure 7]. 

We simulate the creation of a film, or shell on the outside of the tablet as an effect of heating, 
as suggested in |13j . In the model, polymer particles can move towards the surface during thermal 
treatment, resulting in a higher chance of polymer being found on the tablet surface than in the 
interior. The user specifies the thickness of the shell layer (w in the program) and an imbalance 
factor b > 1 that represents the preference of polymer particles to locate in the shell over the 
remainder of the tablet. The distribution of particles is then computed as follows. Let ng denote 
the number of cells (of any type) in the shell layer, and let tit denote the number of cells in the 
remainder of the tablet. To simulate shell formation, a number of np$ — ng miri{&PD, 1} polymer 
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No shell formation Shell formation (b = 2) Complete shell - no release (w = 2) 




FIG. 3: An example of the formation of a polymer shell after thermal treatment. The tablet 
contains 50% polymer. Polymer cells are shown in white, drug and excipient cells are shown in 
black, and gray cells represent mixtures of polymer and other particles. (Left) b = l,w = 1 (default) 
, Center: b = 2, w = 1, giving a shell with some "holes". (Right) b = 2,w = 2, giving a complete 
shell allowing no release. All other parameters are given in Table IT) 



cells are randomly assigned to cells in the shell layer of the tablet. The remaining polymer cells, 
whose number is npx = -Pd(^S + n r) — n PS are randomly assigned to the interior of the tablet. 
All unoccupied tablet cells (including cells in the shell layer) are then reshuffled, and are assigned 
to drug or excipient states, according to their prescribed fractions in the tablet. It should be 
emphasized that setting b = 1 results in no shell formation. See Figure [3] for an illustration of 
simulated shell formations, and Figure [5J top left panel, for a comparison of release curves with and 
without shells. 

Experimental release curves at various polymer concentrations, curing temperatures and com- 
pression pressures show two distinctive characteristics that must be matched by a mathematical 
model. These are the initial slope of the release curve, and the amount of drug that is not released 
after a certain period of time, i.e. the "trapped" drug mass, when the release curve has flattened. 
In order get the observed range of trapped drug mass, the model must reflect the deformation of 
particles under pressure and the fusing of polymer particles under thermal treatment (similar to 
sintering). We emulate these effects by refining the original grid and by allowing polymer particles 
to move into adjacent grid elements, increasing tortuosity and the likelihood of trapping. First, 
each grid element is subdivided into 8 subcells of equal size. If the grid element was in state D 
or X, corresponding to drug or excipient particles, then the states of half of these subcells are 
changed to E, empty. The new empty, state E subcells are located in opposing diagonal elements. 
This subdivision of the cells corresponds to the breaking down of brittle particles under pressure. 
Since polymer particles do not break under pressure, they are subjected to a different "heating" 
rule. Thermal treatment is modeled by swapping void subcells and polymer subcells if a polymer 
particle is located adjacent to a drug or excipient particle. If a polymer particle is adjacent to 
several non-polymer particles, one of these neighbors is chosen randomly and subcells are swapped. 
The processes of subdivision and subcell swapping are shown graphically in Figure [4] The result is 
a virtual tablet that resembles an actual tablet after thermal treatment. 

After initialization, the simulation domain is a cylindrical subregion of a cubic lattice, where each 
grid cell has side length equal to L. If compression is used, then the side length is half the size 
provided by the user L = y/2. The update rules are applied to this computational domain, so that 
a grid cell is one of these cells of length L, and the number of cells in one direction of the cubic 
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FIG. 4: Diagram representing the grid subdivision and "swapping" algorithm used to simulate 
polymer deformation and spreading. Each grid element is subdivided into 8 sub-cells, and empty 
(black) sub-cells are created in drug or excipient grid elements (brown). Polymer sub-cells (blue) 
are swapped with empty (black) sub-cells during the simulation of thermal treatment. 



lattice is n/L. See Figure [2] for a schematic overview of the program flow. 



C. Updating rules 



The cell states are updated according to their present status and the following rules. Table [T] lists 
all user-defined parameters and their default values. 

1. Wetting. For a cell of type E we simulate its filling with water from neighboring cells of type 
W due to capillarity. If a cell of type E has at least one neighbor of type W, then its status 
is changed to W, otherwise no action is taken. 

2. Dissolution. A cell of type D or X that is in contact with cells of type W, can be replaced 
by a W cell. Assume that cell x is in state D, the case that a cell is in state X follows the 
obvious modifications. Denote by nw{x) the number of neighboring water cells of cell x and 
by N\v{x) the set of all neighbors of type W. The dissolution capacity of the water in the 
neighborhood of cell x is 



5=1- 



1 



n w (x) 



J2 <y^- 



y£N w (x) 



Note that S e [0, 1]. We interpret S as the propensity of the particle at position x to dissolve, 
and we implement this by allowing the particle at x to begin to dissolve with probability S. 
In the CA algorithm a random number r is drawn uniformly from [0, 1], and if r < S, the 
cell's status is set to W and the concentration of the substance that just dissolved is set to 1, 
while the other concentration remains at 0. 
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3. Diffusion. Cells of type W have two concentrations associated with them, the concentration 
of dissolved drug, and the concentration of dissolved excipient. These concentrations are 
updated according to the following rule, which is a discretized version of the continuous 
diffusion equation 

c(x, t m+ i) = c(x, t m ) + D \ ^ c(y, t m ) - n w (x)c(x, t m ) 

\yeN w (x) 

where c(x, t m ) denotes the concentration in cell x in the updated tablet and D is the diffusion 
rate of the dissolved substance in cell x relative to the expected time that it takes the fastest 
dissolved substance to diffuse across one grid element. More precisely, 

Di 

D = , 

6 lH(lx{i}^ ru ^, Dexcipient} 

and i stands for "drug" or "excipient". This choice ensures that the new value satisfies 
< c(x) < 1, see Appendix [B] for more details. 

4. Transport. Once dissolved drug or excipient reaches the boundary of the tablet, it must diffuse 
into the surrounding medium. The concentrations in the grid cells surrounding the tablet are 
initially set to zero. As molecules diffuse out from the tablet, the concentrations of drug and 
excipient are updated in the surrounding layer of water cells. Since the tablet is assumed to sit 
in a large amount of fluid, we can assume that concentrations far from the tablet are relatively 
constant. Thus we only model the concentrations in water cells adjacent to the tablet, we call 
the "boundary layer" . In order to simulate different stirring mechanisms, the concentrations 
in cells in the boundary layer are multiplied by a factor f 6 [0, 1] which is provided by the 
user. 

Since particles dissolve and diffuse at different rates, the real time interval corresponding to one 
iteration step, At, is calculated prior to execution of the main loop. Based on rates provided by 
the user, At is set to be the smallest of the four quantities. 

• The average time it takes for a drug respectively excipient particle of mass corresponding to 
one subcell to dissolve. 

• The average time it takes for a dissolved drug respectively excipient particle to diffuse across 
one subcell. 

All other rates are then computed in units of At. Thus, for example, if the average time it takes 
for a particle to dissolve is 2 minutes, and At = 0.6 s, which is the time it takes for a dissolved 
drug molecule to diffuse across one grid element, then the dissolution step takes, on average, ^ = 
200 diffusion steps. We incorporate this variable dissolution rate by decreasing the probability of 
dissolution by a factor of 200, so that, on average, it will take 200 iterations of the algorithm for 
the particle to start to dissolve. See Table [A] for the default values of o and 1, the average increase 
in time to dissolution of drug and excipient particles, respectively. 



D. Release profiles and CA output 



After all cell states have been updated, the simulation clock is advanced it by At, and we calculate 
the percent drug remaining in the tablet. We assume that all of the drug in a D cell can be dissolved 
in a water cell. We will return to this assumption in Section [TV] 
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Let nn(t) be the number of grid cells in state D at time t, so that no(0) is the number of grid 
cells initially in state D. Similarly, let W(f) denote the set of all wet cells in the tablet at time t. 
fraction of the initial drug load in the tablet at time t is 

U ~ n D (0) 

The fraction released at time t is the complement of the fraction of drug in the tablet at that time: 

+E xe w(t) c ( x >*) 



R(t) = 1 - M(t) = 1 - 



n D (0) 



The simulation is stopped once a user-specified time has been reached, or when drug is no longer 
released, i.e. when R'(t) is (approximately) zero. The algorithm detects that drug release has 
stopped when the fraction released has not changed by more than 1CP 5 in 100 consecutive steps. 



III. RESULTS 



The algorithm above has been implemented using C++ in the package celldif f and is available 
at |15j . The simulation time and the fraction of drug released up to that time are exported in ascii 
format and can easily be read using open source languages such as GNUPLOT or SCILAB. In addition, 
the program gives the user the option to export the status of all cells at selected time points as 
an ascii file, allowing for better insight into the wetting and dissolution processes, see Figure [6j 
The state of the simulated tablet suggests explanations for further experimental observations. For 
example, Lemaire et al. [I] asked what is the cause of the differences between the initial release 
of drug and the later release. The CA model describes three phases, namely the uptake of water 
into the tablet, the dissolution of solid drug and excipient, which causes the formation of pores and 
finally the diffusion of dissolved drug out of the tablet. Figure [6] shows an interior slice of a tablet 
8 mm in diameter during the initial 8 minutes of simulation. The water enters the tablet rapidly, 
dissolving the soluble particles of drug and excipient. The interface between solid and dissolved 
drug is approximately circular in cross-section, and with our choice of parameters, the radius of 
the dry region is decreasing at approximately 160 /im min . Such predictions can be tested by 
immersing dry tablets in a dyed solution for different periods of time. 

In Figure [7] we show the simulated release profile of a drug from a matrix tablet of n = 0.8 cm 
diameter with parameters as in Table |lj There is strong qualitative agreement with experimental 
release profiles published in the literature, see e.g. Figure 7], [321 Figures 1-4], (HI Figure 2] and 
[T5I Figures 1-3]. The agreement here is understood as the presence of the following features. 

• Higher polymer fractions in the powder mixture result in partial release, with the released 
fraction decreasing as the polymer fraction increases. 

• The early phase of the release is convex, while the later phase is concave. 

• The release occurs over approximately 8 h. 

We can study the influence of the individual parameters on the shape of the release profiles by 
varying them one at a time, see Figure [8] An increase of the parameter b (Figure [8] top left panel), 
the increased clustering of polymer particles in the outer shell, results in a delayed drug release. In 
the case of higher polymer fractions, this can amount to an almost complete trapping of the drug. 
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FIG. 5: Simulated tablets before and after compression. Left and middle: a slice of the tablet 
before and after compression. Polymer cells are lighter, drug and excipient are gray, and boundary, 
wet or void cells are black. The lightest cells (middle panel) are cells in which polymer and drug 
are mixed. Right: Release curves for a tablet with compression and without compression. Each 
tablet contains 50% polymer. Other parameters are as given in Table [TJ 



The stirring of the surrounding fluid [8] decreases the concentration of dissolved drug outside the 
tablet and thereby facilitates the diffusion, although this effect is of somewhat limited size. Clearly, 
this is of importance mostly for experiments in rotating disk apparatuses and similar devices. An 
increase of either the diffusion rate k or the dissolution rate, o results in a faster release of the drug. 
We remark here that the dissolution parameters o and 1 affect the rate of erosion of the tablet, and 
hence their effect is most noticeable at the initiation of release. The parameter k affects the rate 
at which drug diffuses from the tablet, and so has a more pronounced effect over the entire release 
curve. 



IV. DISCUSSION 



From our simulations we gain insight on the role that the individual parameters play and how 
they influence the drug release profiles. Of course, it has to be borne in mind that some parameters 
cannot be manipulated in the process of tablet formulation and fabrication, such as dissolution and 
diffusion rates. Roughly speaking, we can group the parameters into two categories, namely those 
that determine the size of the tablet and those that determine its topology. In the former category 
are the diffusion rate u, the dissolution rate o and, naturally, the tablet diameter n. In the latter 
category are the width of the melted polymer shell w and the imbalance factor b that models how 
much more the polymer particles are clustered in the melting zone at the edge of the tablet. 

A natural question is the powder composition necessary to produce an "ideal" tablet, with a 
largely linear release profile and a almost complete release of the (expensive) drug. According 
to our simulation results and consistent with experimental evidence |I3j . thermal treatment after 
compression causes drug to be trapped, except at very low polymer fractions. At these low polymer 
fractions, however, drug is released very quickly since most of the tablet erodes. Therefore, in order 
to get complete release of the drug over 8 hours, with a nearly linear release profile, the model 
suggests using a polymer fraction between 40 and 45%, with as little deformation and fusing of 
the polymer particles as possible. We illustrate this in Figure |9j which shows simulations in which 
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FIG. 7: Simulation of drug release from a tablet of n = 0.8 cm diameter using the parameters in 

Table [TJ 



the compression sub-routine is disabled (flag d is set to 0). Note that in this model we assume 
brittle drug and excipient, such as indomethacin and lactose, and a deformable, non-soluble, non- 
permeable polymer. If the properties of the particles are changed, then the model must be modified 
accordingly. For a release profile close to the ideal it is desirable to put a shell of polymer on the 
tablet without completely enclosing the interior. In our model, we set b = 2, but we do not apply 
compression or heating. With 30% polymer, we get a close to ideal release over 8 hours, and higher 
polymer concentrations extend the release profile over longer time intervals. See Figure [9] (right 
panel) . 

Of the possible drug release mechanisms, our model mimics the transport of drug molecules 
through water filled pores. As time passes, degradation of the excipient and drug particles occurs 
and many crevices and channels are formed. The diffusion through these pores is highly dependent 
on the polymer structure and consequently also dependent on the processes that alter pore formation 
and closure pQ . As it can be seen in Figure [5] the spreading and fusing of polymer particles after 
compression and thermal treatment changes the tablet structure and modifies the tortuosity, pore 
formation and the resulting release curves. This structural change is crucial to reproduce the 
trapping of drug at high polymer concentrations observed in experimental data. Without it, the 
release obtained at high polymer concentrations did not match that of the experimental data. 
The initial burst release observed in Figure [7] for all polymer densities, can be attributed to drug 
particles easily accessible by hydration on the surface [HHZj- Thus, to match the very low rate 
of drug release present at high curing temperatures and high polymer content, we implement the 
formation of a polymer shell on the surface of the tablet. As described in [T3], the effect of thermal 
treating causes a better coalescence of polymer particles with a decreased porosity and a smoother 
surface. By altering the morphology of the outer layer, we successfully control the permeability and 
eliminate the initial burst, providing a more faithful picture of the experimental results. 

To include other release mechanisms such as diffusion through the polymer or degradation of 
the polymer in our mathematical model, it will be necessary to describe the diffusion coefficient as 
a function of other parameters. The diffusion rate through a polymer is highly dependent on its 
physical state. For some polymers, there is a glass temperature transition in which the polymer 
changes from brittle to rubber-like. In this case, the diffusion rate can change by several orders 
of magnitude pQ. For biodegradable polymers, modeling the dissolution of the polymer and the 
corresponding matrix degradation might require a non-constant diffusion coefficient that varies in 
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FIG. 9: (Left) Simulated release curves for tablets containing 40, 45 and 50 % polymer, modeling 
no compression or thermal treatment. The compression flag, d, is set to 0. All other parameters 
are in Table |TJ (Right) Simulated release curves for tablets containing 30, 40 and 45 % polymer, 
with no compression (d = 0), with an imposed shell of polymer, b = 2. Comparison is made to a 
tablet with 40% polymer and no shell. The release curves are close to the "ideal", depending on 

the desired release interval. 
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proportion to the local fraction of polymer in the tablet [T8] , 

The role of the excipient in our model is merely that of a space filler, which is not to say that it is 
not important. Dissolving excipient opens new channels through which water can reach the drug, 
and through which the drug can diffuse out of the tablet. It should be stressed that our assumption 
that drug and excipient dissolve independently of each other in water is a simplification. A water cell 
whose water is saturated with respect to excipient should have a very small capacity to dissolve drug 
and vice versa. Furthermore, the current assumption that the entire mass of drug in one grid cell 
can be dissolved in the amount of water in a grid cell of the same volume is another simplification. 
A more realistic implementation of the model would incorporate a gradual dissolution process, so 
that only a fraction of the drug in a D cell is dissolved in one simulation step. 

One simulation representing 8 hours with the baseline parameters given in Table [I] typically 
takes one hour on a MacPro workstation (6-Core Intel Xeon at 2.93 GHz). The algorithm is easily 
parallelizable, which is one of the positive features of a CA model. This is work in progress. 



A sample call of the code is as follows: user-supplied arguments are preceded by a hyphen, and 
successive arguments are separated by spaces. For example, the call 

./celldiff -n 0.8 -p 0.3 -c30000 -tlOO -sTabletData.txt 

calls the program simulating a tablet with diameter n = 0.8 cm, made up of 30% polymer and 
10% drug (the default). It simulates the evolution of the tablet over 3 • 10 4 s and exports the full 
state data every 100 time steps, storing it in the text file TabletData.txt. Table |T] lists the model 
parameters and their default values in the program, celldiff. 



Let c(x,t) be the concentration of a substance at a point x € R 3 at time t € K. The diffusion 
equation in three dimensions is 



where A denotes the usual Laplace operator in three dimensions and D is the diffusion constant 
(with units Length 2 /Time). After discretization of space on a cubic lattice with lattice constant L 
and of time, so that x — (i,j,k) and t m — mAt, the discretized version of the diffusion equation 
becomes 



where N(x) is the von Neumann neighborhood of x. In the CA model, we choose At to be the 
expected time it would take for a dissolved particle to diffuse across one grid element, of side length 
L. Since the variance of the diffusion process, X(t), is given by i?[A 2 ] = 6Dt, we choose the 
time step, At so that I? = 6DAt and hence At = gjj, where D = max{Z?,j rtig , D excipient } is the 



Appendix A: Sample call and parameter values 



Appendix B: Derivation of the diffusion update rule 




(Bl) 
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Name 


User Input 


Meaning 


Default Value 


Units/Range 


diameter 


n 


diameter of the computa- 
tional domain 


0.8 


cm 


maxtime 


c 


maximum time to simulate 
before halting 


30000 


s 


polymer 


P 


polymer concentration 


0.3 


[0,1] 


drug 


g 


drug concentration 


0.1 


[0,1] 


tablet height 


h 


height of tablet as a ratio 
of diameter 


0.23 


[0, i] 


release file 


r 


filename for exporting re- 
lease curve 


time stamped name 




state file 


s 


filename for exporting full 
state data 


time stamped name 




stateperiod 


t 


how often to export full 
state data 


(never) 




compress 


d 


flag to indicate whether 
or not compression and 
thermal treatment is 
implemented 


1 


1 or 


seed 


e 


random number seed 


47 


N 


asciiperiod 


a 


frame period in iterations 
for drawing ascii animation 


1 


N 


drugdissprob 





dissolution probability 
scale for drug 


10 _a 




exdissprob 


1 


dissolution probability 
scale for excipient 


10" 3 




polyshellbalance 


b 


distribution of polymer in 
+Vip qTipII 

LUC ollLll 


1 


> 1 


rpmnvfil vatp 


f 


rpmnvfll factor for pnnppn- 

trations in boundary cells 


0.95 


fo 11 


nographics 


X 


set > to run in non- 
interactive, text-only mode 


1 




drugdiffrate 


u 


physical rate of drug 
diffusion 


7- 10" 6 


cm 
1? 


exdiffrate 


k 


physical rate of excipient 
diffusion 


7- 10" 6 


cm 

1? 


cellsize 


y 


physical size of a cell 


0.01 


cm 


polyshellwidth 


w 


width of the polymer shell 
(has no effect when 6=1, 
the default) 


1 


cell size 



TABLE I: Default values of the parameters in celldif f . 



maximum of the diffusion rates of the drug and excipient. Substituting this value of At into the 



discretized diffusion equation, equation (Bl|, and rearranging gives 



c(x, t m+ i) = c(x, t m ) + D [ ^ c(y, t m ) - 6c(x, t % 



(B2) 
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B DERIVATION OF THE DIFFUSION UPDATE RULE 



A 

where D = jr— — , where i is "drug" or "excipient" , depending on the concen- 

D IiltXX\U d rU g , U excipient J 

tration that is updated. In this model, certain neighbors might contain solid particles of drug or 
excipient, or polymer. We impose zero flux conditions on the boundaries to these solid neighbors, 
so that only terms involving neighboring cells in the "wet" (W) state are counted in the expression 
in equation (B2). This gives the diffusion update rule described above. 
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